Let’s revise the priors to be more consistent with current practice. Specifically,
rm(list = ls())
gc()
modelName <- "atorv1"
scriptName <- paste(modelName, "Rmd", sep = ".")
fitModel <- FALSE
## Relative paths assuming the working directory is the script directory
## containing this script
scriptDir <- getwd()
projectDir <- dirname(scriptDir)
dataDir <- file.path(projectDir, "data")
modelDir <- file.path(projectDir, "model")
outDir <- file.path(modelDir, modelName)
figDir <- file.path(projectDir, "deliv", "figure", modelName)
tabDir <- file.path(projectDir, "deliv", "table", modelName)
invisible(dir.create(figDir, recursive = TRUE))## Warning in dir.create(figDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/figure/atorv1' already exists
invisible(dir.create(tabDir, recursive = TRUE))## Warning in dir.create(tabDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/table/atorv1' already exists
.libPaths("lib")
library(metrumrg)## Loading required package: reshape
## Loading required package: lattice
## Loading required package: XML
## Loading required package: MASS
## Loading required package: grid
## metrumrg 5.57
## enter "?metrumrg" for help
suppressMessages(library(rstan))
suppressMessages(library(bayesplot))
suppressMessages(library(tidyverse))
library(gridExtra)##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(parallel)
library(knitr)
library(kableExtra)##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(qapply)## Loading required package: rlecuyer
##
## Attaching package: 'qapply'
## The following object is masked from 'package:metrumrg':
##
## qstat
library(PKPDmisc)
knitr::opts_chunk$set(
comment = '.',
fig.height = 5,
fig.width = 9,
eval.after = 'fig.cap',
message = FALSE,
fig.path = file.path(figDir, paste(modelName, "_", sep = ""))
)
## Go back to default ggplot2 theme that was overridden by bayesplot
theme_set(theme_gray())
source(file.path(scriptDir, "tools", "stanTools2.R"))
source(file.path(scriptDir, "tools", "functions.R"))
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(11191951)
set.seed(10271998) modelFile <- file.path(modelDir, paste(modelName, "stub.ctl", sep = ""))
modelText <- readLines(modelFile)
cat(modelText, sep = "\n"). $PROBLEM RUN# atorv1
. $INPUT C ID TIME AMT DV CMT DOSE EVID MDV ADDL II SS FORM WT PER BLQ
. ;$INPUT C ID SID=DROP STDY=DROP TIME DV DVR=DROP DVC=DROP AMT AMTO=DROP CMT TYPE=DROP
. ; DOSE DOSN=DROP NTIM=DROP EVID MDV ADDL II SS FORM WT SEX=DROP AGE
. ; HT=DROP RACE=DROP PER BLQ TRCD=DROP ANLT=DROP TXT=DROP
. ;$DATA ../../../data/atorv1172atvb.csv IGNORE=(C='C', BLQ.EQ.1)
. $DATA ../../../data/atorvWrkShop.csv IGNORE=(C='C', BLQ.EQ.1)
.
. $ABBR DECLARE INTEGER FIRST_WRITE_PAR
. $ABBR DECLARE INTEGER FIRST_WRITE_IPAR
.
. $SUBROUTINES ADVAN4 TRANS4 OTHER = ../../extrasend.f90
.
. $PK
. ; Request extra information for Bayesian analysis.
. ; An extra call will then be made for accepted samples
. include '/opt/NONMEM/nm74/util/nonmem_reserved_general'
. BAYES_EXTRA_REQUEST=1
.
. nThin = THETA(8)
.
. VWT = LOG(WT/70) ; normalized to 70 kg adult
. MU_1 = THETA(1) + VWT * 0.75 ; CL
. MU_2 = THETA(2) + VWT ; V2
. MU_3 = THETA(3) + VWT * 0.75 ; Q
. MU_4 = THETA(4) + VWT ; V3
. MU_5 = THETA(5) ; ka
. MU_6 = THETA(6) ; F1chew
. MU_7 = THETA(7) ; SD
.
. " CALL EXTRASEND()
.
. ;ATORVASTATIN
.
. CL = EXP(MU_1 + ETA(1)) ; atv cl
. V2 = EXP(MU_2 + ETA(2))
. Q = EXP(MU_3 + ETA(3))
. V3 = EXP(MU_4 + ETA(4))
. ;deltaKA = EXP(MU_5 + ETA(5))
. KA = EXP(MU_5 + ETA(5))
. F1chew = EXP(MU_6 + ETA(6))
. SD = EXP(MU_7 + ETA(7))
.
. ; Constrain KA > lambda1 (typical value) to enhance identifiability
. ;T1 = CL/V2
. ;T23 = Q/V2
. ;T32 = Q/V3
. ;L1 = ((T1+T23+T32)+SQRT((T1+T23+T32)**2-4*T1*T32))/2
. ;L2 = ((T1+T23+T32)-SQRT((T1+T23+T32)**2-4*T1*T32))/2
.
. ;KA = deltaKA + L2
. S2 = V2/1000 ; DOSE IN uM & CONC IN nM/L
. F1 = 1.0 ; TABLET AS REFERENCE
. IF(FORM.EQ.2) F1 = F1chew ; CHEWABLE F1
.
. IF(BAYES_EXTRA==1 .AND. NIREC==1 .AND. NDREC==1 .AND. &
. ITER_REPORT>0 .AND. &
. ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN
. IF(FIRST_WRITE_PAR==0) THEN
. " OPEN(unit=52,FILE='./par.txt')
. FIRST_WRITE_PAR=1
. ENDIF
. " WRITE(52,'(I12,1X,20(1X,1PG19.10E3))') ITER_REPORT, &
. " OMEGA(1,1), OMEGA(2,1), OMEGA(2,2), &
. " THETA(1), THETA(2), THETA(3), THETA(4), THETA(5), THETA(6), THETA(7)
. ENDIF
. IF(BAYES_EXTRA==1 .AND. NDREC==1 .AND. ITER_REPORT>0 .AND. &
. ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN
. IF(FIRST_WRITE_IPAR==0) THEN
. " OPEN(unit=50,FILE='./ipar'//TFI(PNM_NODE_NUMBER)//'.txt')
. FIRST_WRITE_IPAR=1
. ENDIF
. " WRITE(50,'(I12,1X,F14.0,10(1X,1PG19.10E3))') ITER_REPORT, ID, &
. " ETA(1), ETA(2)
. ENDIF
.
.
. $ERROR
.
. LOQ = 0.45 ; nM/L
. IPRED = F
. DUM = (LOQ-IPRED)/(SD*IPRED)
. CUMD=PHI(DUM)
. IF(BLQ.EQ.0)THEN
. F_FLAG=0
. Y = F*(1+SD*ERR(1)) ; ATV
. ELSE
. F_FLAG=1
. Y=CUMD
. ENDIF
.
. ; Code below into chains
## create initial estimates
geninit <- function(){
paste(c(
"; Initial values of THETA",
"$THETA",
paste(rnorm(1, log(1000), 0.5), " ; TOTAL CL ATV = THETA(1)"),
paste(rnorm(1, log(600), 0.5), " ; V2 = THETA(2)"),
paste(rnorm(1, log(130), 0.5), " ; Q = THETA(3)"),
paste(rnorm(1, log(1500), 0.5), " ; V3 = THETA(4) "),
paste(rnorm(1, log(1.2), 0.5), " ; KA= THETA(5) "),
paste(rnorm(1, log(1), 0.25), " ; THETA(6) relative F1 CHEWABLE = THETA(5) "),
paste(log(0.2), " ; THETA(7) - PRO RES ERROR"),
"$OMEGA BLOCK(2) ;INITIAL values of OMEGAs",
"0.1003 ; cl",
".06034 0.8477 ; v2",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$SIGMA ;Initial value of SIGMA",
"1 FIX; pro"),
sep = "\n")
}
## Set parameters of the prior distribution
priors <- paste(c("$PRIOR NWPRI",
"$THETAP ; Prior information of THETAS",
paste("(", 5.65, " FIX) ; THETA(1) TOTAL CL ATV"),
paste("(", 5.60, " FIX) ; THETA(2) V2"),
paste("(", 4.29, " FIX) ; THETA(3) Q"),
paste("(", 7.07, " FIX) ; THETA(4) V3"),
paste("(", -1.33, " FIX) ; THETA(5) KA"),
paste("(", log(1), " FIX) ; THETA(6) F1 chewable"),
paste("(", log(0.2), " FIX) ; THETA(7) log(sigma)"),
"$THETAPV BLOCK(7) ; variances for priors on THETAS (var-cov)",
"2 FIX ; CL weakly informative",
"0.00 2 ; V2 weakly informative",
"0.00 0.00 2 ; Q weakly informative",
"0.00 0.00 0.00 .0156 ; V3 SE**2 informative",
"0.00 0.00 0.00 0.00 .00684 ; KA SE**2 informative",
"0.00 0.00 0.00 0.00 0.00 1 ; F1 chewable weakly informative",
"0.00 0.00 0.00 0.00 0.00 0.00 1 ; log(sigma) weakly informative",
## Parameters LKJ / lognormal prior; off-diagonal and OMEGAP ignored
"$OMEGAP BLOCK(2) ; prior information for OMEGA",
".1003 FIXED ; cl",
".06034 .8477 ; v2" ,
"$OMEGAPD (4.0, FIXED) ; df for OMEGA prior"
## Parameters for inverse Wishart prior
## "$OMEGAP BLOCK(2) ; prior information for OMEGA",
## ".1003 FIXED ; cl",
## ".06034 .8477 ; v2" ,
## "$OMEGAPD (4.0, FIXED) ; df for OMEGA prior"
),
sep = "\n")
## Table statements
tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE WT",
"NOPRINT ONEHEADER FILE=./.tab",
"$TABLE ID WT CL V2 Q KA V3 ETA1 ETA2 PER FORM",
"NOPRINT ONEHEADER FILE=./par.tab"),
sep = "\n")
##tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE AGE WT",
## "NOPRINT ONEHEADER FILE=./.tab",
## "$TABLE ID WT CL V2 Q KA V3 ETA1 ETA2 PER FORM",
## "NOPRINT ONEHEADER FILE=./par.tab"),
## sep = "\n")nChains <- 4
nPost <- 250 ## Number of post-burn-in samples per chain after thinning
nBurn <- 250 ## Number of burn-in samples per chain after thinning
nThin <- 1
seed = sample(10000, size = nChains)
if(fitModel){
if(!file.exists(file.path(modelDir, modelName)))
dir.create(file.path(modelDir, modelName))
runChains <- mclapply(1:nChains, runChain,
modelName = modelName, modelDir = modelDir,
priors = priors, tables = tables,
nPost = nPost, nBurn = nBurn, nThin = nThin,
seed = seed, print = 100,
## LKJ prior for IIV correlation matrix with oarameter = OLKJDF
OLKJDF = 1,
## lognormal prior for IIV SDs
## log(sqrt(Omega(i))) ~ Normal(log(sqrt(OmegaPrior(i))), 1/OVARF)
OVARF = 1, AUTO = 2,
NUTS_DELTA = 0.8,
grid = FALSE,
method = "NUTS",
## pe = "orte 16",
mode = "nonpara")
# NONR(run = paste(modelName, ".", 1, sep = ""),
# command = "/opt/NONMEM/nm74/nmqual/autolog.pl",
# project = file.path(modelDir, modelName),
# grid = FALSE,
# wait = FALSE,
# diag = FALSE,
# fdata = FALSE,
# purge = FALSE,
# checkrunno = TRUE)
}dataFile <- "atorvWrkShop.csv"
nmData <- read.csv(file.path(dataDir, dataFile), as.is = TRUE)chains <- 1:nChains
##chains <- c(1, 4)
nChains2 <- length(chains)
## Read in posterior distributions by chain and add a column called chain
popParameters <- map(chains, function(thisChain){
param <- data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""),
"par.txt"), data.table = FALSE)
names(param) <- c("iteration", "OMEGA11", "OMEGA21", "OMEGA22",
"THETA1", "THETA2", "THETA3", "THETA4", "THETA5",
"THETA6", "THETA7")
param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
mutate(sample = 1:n())
## Calculate more interpretable paraneters
popParameters2 <- popParameters %>%
mutate(CLHat = exp(THETA1),
V2Hat = exp(THETA2),
QHat = exp(THETA3),
V3Hat = exp(THETA4),
## k10 = CLHat / V2Hat,
## k12 = QHat / V2Hat,
## k21 = QHat / V3Hat,
## lambda = ((k10 + k12 + k21) -
## sqrt((k10 + k12 + k21)^2 - 4 * k10 * k21)) / 2,
## KAHat = exp(THETA4) + lambda,
## KAHat = exp(THETA4 + lambda),
KAHat = exp(THETA5),
F1Chewable = exp(THETA6),
sigma1 = exp(THETA7)) %>%
select(starts_with("OMEGA"),
CLHat, V2Hat, QHat, V3Hat, KAHat, F1Chewable, sigma1, chain, iteration, sample)
## Convert to 3-D array with dims = {iterations, chains, parameters}.
popParArray <- array(double(nPost * nChains2 * (ncol(popParameters2) - 3)),
dim = c(nPost, nChains2, ncol(popParameters2) - 3),
dimnames = list(NULL, NULL, setdiff(names(popParameters2), c("chain", "iteration", "sample"))))
for(iChain in 1:nChains2){
popParArray[,iChain,] <- popParameters2 %>%
filter(chain == chains[iChain],
iteration > 0) %>%
select(-chain, -iteration, -sample) %>%
as.matrix
}indParameters = map(chains, function(thisChain){
files <- list.files(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = "")))
iparFiles <- files[grepl("^ipar", files)]
param <- map(iparFiles, function(thisFile){
data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""),
thisFile), data.table = FALSE)
}) %>%
bind_rows
names(param) <- c("iteration", "ID", "ETA1", "ETA2")
param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
arrange(chain, iteration, ID) %>%
left_join(popParameters)options(bayesplot.base_size = 12,
bayesplot.base_family = "sans")
color_scheme_set(scheme = "brightblue")
myTheme <- theme(text = element_text(size = 12),
axis.text = element_text(size = 12))
plot_mcmcHistory <- mcmcHistory(popParArray,
nParPerPage = 5, myTheme = myTheme)
plot_mcmcDensityByChain <- mcmcDensity(popParArray,
nParPerPage = 16, byChain = TRUE,
myTheme = theme(text = element_text(size = 12),
axis.text = element_text(size = 10)))
plot_mcmcDensity <- mcmcDensity(popParArray, nParPerPage = 16,
myTheme = theme(text = element_text(size = 12),
axis.text = element_text(size = 10)))
plot_mcmcHistory. [[1]]
Pediatric atorvastatin PK example: MCMC trace plots.
.
. [[2]]
Pediatric atorvastatin PK example: MCMC trace plots.
plot_mcmcDensityByChain. [[1]]
Pediatric atorvastatin PK example: Univariate marginal densities by chain.
plot_mcmcDensity. [[1]]
Pediatric atorvastatin PK example: Univariate marginal densities.
mcmc_pairs(popParArray[,,!grepl("^OMEGA", dimnames(popParArray)[[3]])])Pediatric atorvastatin PK example: Bivariate marginal densities.
mcmc_pairs(popParArray[,,grepl("^OMEGA", dimnames(popParArray)[[3]])])Pediatric atorvastatin PK example: Bivariate marginal densities.
caption <- paste("Pediatric atorvastatin PK example:",
c(rep("MCMC trace plots.", length(plot_mcmcHistory)),
rep("Univariate marginal densities by chain.",
length(plot_mcmcDensityByChain)),
rep("Univariate marginal densities.",
length(plot_mcmcDensityByChain)),
rep("Bivariate marginal densities.", 2)))
ptable <- monitor(popParArray,
warmup = 0, print = FALSE) %>%
as.matrix %>%
formatC(3) %>%
as.data.frame
write.csv(ptable, file = file.path(tabDir, paste(modelName,
"ParameterTable.csv", sep = "")))
ptable %>%
rename(SEmean = se_mean, SD = sd, pct2.5 = "2.5%", pct25 = "25%", median = "50%",
pct75 = "75%", pct97.5 = "97.5%", Neff = "n_eff") %>%
mutate(parameter = rownames(.), "95% CI" = paste("(", pct2.5, ", ", pct97.5, ")",
sep = "")) %>%
select(parameter, mean, SD, median, "95% CI", Neff, Rhat) %>%
kable(caption = "Summary of model parameter estimates.") %>%
kable_styling(bootstrap_options = "striped", full_width = F)| parameter | mean | SD | median | 95% CI | Neff | Rhat |
|---|---|---|---|---|---|---|
| OMEGA11 | 0.416 | 0.171 | 0.419 | (0.195, 0.769) | 3 | 1.56 |
| OMEGA21 | 0.0898 | 0.128 | 0.125 | (-0.213, 0.305) | 37 | 1.11 |
| OMEGA22 | 294 | 508 | 1.04 | (0.48, 1.17e+03) | 2 | 1.53 |
| CLHat | 590 | 151 | 618 | ( 373, 842) | 3 | 1.57 |
| V2Hat | 2.44e+03 | 2.36e+03 | 1.17e+03 | ( 667, 6.52e+03) | 2 | 1.53 |
| QHat | 394 | 213 | 293 | ( 188, 755) | 2 | 1.52 |
| V3Hat | 1.6e+03 | 200 | 1.6e+03 | (1.37e+03, 2.01e+03) | 4 | 1.33 |
| KAHat | 0.27 | 0.0396 | 0.281 | (0.209, 0.329) | 2 | 1.53 |
| F1Chewable | 0.896 | 0.179 | 0.87 | (0.675, 1.21) | 374 | 1.55 |
| sigma1 | 0.258 | 0.0494 | 0.235 | (0.208, 0.341) | 2 | 1.53 |
## For now just use the NONMEM data set
design <- nmData %>%
filter(C == ".")library(mrgsolve)
indCode =
'
$PARAM
THETA1 = 5.65
THETA2 = 5.60
THETA3 = 4.29
THETA4 = -1.33
THETA5 = 7.07
THETA6 = 1
THETA7 = -3.2189
WT = 70
FORM = 1
ETA1 = 0
ETA2 = 0
$MAIN
double VWT = log(WT/70);
double CL = exp(THETA1 + (VWT * 0.75) + ETA1 + ETA(1));
double V2 = exp(THETA2 + VWT + ETA2 + ETA(2));
double Q = exp(THETA3 + (VWT * 0.75));
double V3 = exp(THETA5 + VWT);
//double T1 = CL/V2;
//double T23 = Q/V2;
//double T32 = Q/V3;
//double L2 = ((T1+T23+T32)-sqrt(pow(T1+T23+T32, 2)-4*T1*T32))/2;
double KA = exp(THETA4);
F_GUT = 1.0; // TABLET AS REFERENCE
if(FORM == 2) F_GUT = exp(THETA6); // CHEWABLE F1
$PKMODEL cmt = "GUT CENT PERIPH", depot = TRUE
$SIGMA 1
$OMEGA @block
0 0 0
$TABLE
double cPredHat = 1000 * (CENT / V2);
double cPred = cPredHat * (1 + exp(THETA7) * EPS(1));
$CAPTURE cPred
'
## Compile model
atvMod <- mcode("ATVmodel", indCode)
if(FALSE){
## Test runs
out <- atvMod %>%
data_set(design) %>%
carry_out(DOSE) %>%
mrgsim
plot(out, cPred ~ TIME | factor(DOSE), scales = "same")
out <- atvMod %>%
data_set(design) %>%
param(popParameters2[1,] %>%
rename(CL = CLHat,
V2 = V2Hat,
Q = QHat,
V3 = V3Hat,
KA = KAHat,
F1 = F1Chewable)) %>%
carry_out(DOSE) %>%
mrgsim
plot(out, cPred ~ TIME | factor(DOSE), scales = "same")
}samples <- unique(indParameters$sample)
samples <- 1:8
indSim1 <- mclapply(samples,
function(thisSample){
atvMod %>%
data_set(design) %>%
idata_set(indParameters %>%
filter(sample == thisSample)) %>%
carry_out(sample, ID, TIME, DV, EVID, DOSE, FORM) %>%
mrgsim %>%
as.tbl
}) %>%
bind_rowssamples <- unique(popParameters$sample)
popSim1 <- mclapply(samples,
function(thisSample){
atvMod %>%
data_set(design) %>%
param(popParameters %>%
filter(sample == thisSample)) %>%
omat(matrix(with(popParameters %>%
filter(sample == thisSample),
c(OMEGA11, OMEGA21, OMEGA21, OMEGA22)), ncol = 2)) %>%
carry_out(sample, ID, TIME, DV, EVID, DOSE, FORM) %>%
mrgsim %>%
as.tbl
}) %>%
bind_rows## Combined plots
indPred <- indSim1 %>%
group_by(ID, TIME, EVID, DOSE, FORM) %>%
summarize(lbInd = quantile(cPred, probs = 0.05, na.rm = TRUE),
medianInd = quantile(cPred, probs = 0.5, na.rm = TRUE),
ubInd = quantile(cPred, probs = 0.95, na.rm = TRUE))
popPred <- popSim1 %>%
group_by(ID, TIME, EVID, DOSE, FORM) %>%
summarize(lbPop = quantile(cPred, probs = 0.05, na.rm = TRUE),
medianPop = quantile(cPred, probs = 0.5, na.rm = TRUE),
ubPop = quantile(cPred, probs = 0.95, na.rm = TRUE))
predAll <- indPred %>%
left_join(popPred) %>%
left_join(design) %>%
mutate(DV = suppressWarnings(as.numeric(DV)))
forms <- sort(unique(predAll$FORM))
plot_PPCPK <- lapply(forms,
function(thisForm){
p1 <- ggplot(predAll %>% filter(FORM == thisForm),
aes(x = TIME, y = DV))
p1 <- p1 +
geom_line(aes(x = TIME, y = medianPop,
color = "population")) +
geom_ribbon(aes(ymin = lbPop, ymax = ubPop,
fill = "population"),
alpha = 0.25) +
geom_line(aes(x = TIME, y = medianInd,
color = "individual")) +
geom_ribbon(aes(ymin = lbInd, ymax = ubInd,
fill = "individual"),
alpha = 0.25) +
scale_color_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1") +
scale_fill_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1")
p1 + geom_point() +
labs(title = paste("formulation =", c("tablet", "chewable")[thisForm]),
x = "time (h)",
y = "atorvastatin plasma concentration (nM)") +
theme(text = element_text(size = 12),
axis.text = element_text(size = 12),
## legend.position = c(0.8, 0.25),
strip.text = element_text(size = 8)) +
facet_wrap(~ ID, scales = "free_y")
})
plot_PPCPK[[1]]
. Warning: Removed 48 rows containing missing values (geom_path).
. Warning: Removed 367 rows containing missing values (geom_point).
[[2]]
. Warning: Removed 227 rows containing missing values (geom_point).
sessionInfo(). R version 3.5.1 (2018-07-02)
. Platform: x86_64-pc-linux-gnu (64-bit)
. Running under: Ubuntu 14.04.5 LTS
.
. Matrix products: default
. BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
. LAPACK: /usr/lib/lapack/liblapack.so.3.0
.
. locale:
. [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
. [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
. [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
. [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
. [9] LC_ADDRESS=C LC_TELEPHONE=C
. [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
.
. attached base packages:
. [1] parallel grid stats graphics grDevices utils datasets
. [8] methods base
.
. other attached packages:
. [1] mrgsolve_0.9.2 PKPDmisc_2.1.1 qapply_1.39
. [4] rlecuyer_0.3-4 kableExtra_1.1.0 knitr_1.25
. [7] gridExtra_2.3 forcats_0.4.0 stringr_1.4.0
. [10] dplyr_0.8.3 purrr_0.3.2 readr_1.3.1
. [13] tidyr_1.0.0 tibble_2.1.3 tidyverse_1.2.1
. [16] bayesplot_1.7.0 rstan_2.19.2 ggplot2_3.2.1
. [19] StanHeaders_2.19.0 metrumrg_5.57 MASS_7.3-51.1
. [22] XML_3.98-1.20 lattice_0.20-38 reshape_0.8.8
.
. loaded via a namespace (and not attached):
. [1] httr_1.4.1 jsonlite_1.6
. [3] viridisLite_0.3.0 modelr_0.1.5
. [5] assertthat_0.2.1 highr_0.8
. [7] stats4_3.5.1 cellranger_1.1.0
. [9] yaml_2.2.0 pillar_1.4.2
. [11] backports_1.1.5 glue_1.3.1
. [13] digest_0.6.21 RColorBrewer_1.1-2
. [15] rvest_0.3.4 colorspace_1.4-1
. [17] htmltools_0.4.0 plyr_1.8.4
. [19] pkgconfig_2.0.3 broom_0.5.2
. [21] haven_2.1.1 scales_1.0.0
. [23] webshot_0.5.1 processx_3.4.1
. [25] generics_0.0.2 withr_2.1.2
. [27] lazyeval_0.2.2 cli_1.1.0
. [29] magrittr_1.5 crayon_1.3.4
. [31] readxl_1.3.1 evaluate_0.14
. [33] ps_1.3.0 nlme_3.1-137
. [35] RcppArmadillo_0.9.800.1.0 xml2_1.2.2
. [37] pkgbuild_1.0.6 data.table_1.12.4
. [39] tools_3.5.1 loo_2.1.0
. [41] prettyunits_1.0.2 hms_0.5.1
. [43] lifecycle_0.1.0 matrixStats_0.55.0
. [45] munsell_0.5.0 callr_3.3.2
. [47] compiler_3.5.1 rlang_0.4.0
. [49] ggridges_0.5.1 rstudioapi_0.10
. [51] labeling_0.3 rmarkdown_1.16
. [53] gtable_0.3.0 inline_0.3.15
. [55] reshape2_1.4.3 R6_2.4.0
. [57] lubridate_1.7.4 zeallot_0.1.0
. [59] stringi_1.4.3 Rcpp_1.0.2
. [61] vctrs_0.2.0 tidyselect_0.2.5
. [63] xfun_0.10